%% This Dynare file solves the baseline model in "Financial Constraints, 
%%  Sectoral Heterogeneity, and the Cyclicality of Investment" by Cooper Howes

%% Add the path that includes Dynare files
% addpath c:\dynare\5.2\matlab

%% Change working directory to folder containing replication materials
% cd 

%% Declare variables

var

%% HH variables

% Borrowers
c_n_b c_d_b d_b h_b psi lambda_b b md_b

% Savers
c_n_s c_d_s d_s h_s lambda md_s

% Aggregates
k d c c_d y mpk_d mpk_n

%% Firm variables
mc_d mc_n mk_d mk_n y_d y_n h_d h_n k_d k_n i_d i_n mu

%% Prices
r w pi p_d p_n pi_d pi_n

%% Exogenous stuff
a money xi nu delta sigma;

% Exogenous shocks
varexo e_m e_a e_nu e_sigma e_xi e_delta;

%% Declare parameters
parameters beta beta_b sigma_bar nu_bar chi phi_n phi_d delta_bar delta_d alpha_d alpha_n xibar rho rho_m omega
phi_pi theta_d theta_n theta_c eta_d eta_n m b_d b_n phi_y hab wag_rid w_markup epsilon_d epsilon_n nu_b;

%% Provide parameter values
beta = 0.99; % Discount rate for savers
beta_b = 0.98; % Discount rate for borrowers
sigma_bar = 1; % CRRA parameter
nu_bar = 4; % Labor disutility scale (savers)
nu_b = 4; % Labor disutility scale (borrowers)
chi = 1; % Labor elasticity
theta_d = 2; % Capital investment adjustment costs for durables sector
theta_n = 2; % Capital investment adjustment costs for nondurables sector
theta_c = 0; % Durable investment adjustment costs for consumers
delta_bar = 0.03; % Depreciation rate for capital
delta_d = 0.02; % Depreciation rate for consumer durables
alpha_d = 1/3; % Capital share
alpha_n = 1/3; % Nondurable capital share
xibar = 0.1; % Steady state borrowing constraint
rho = 0.9; % Persistence parameter
rho_m = 0; % Monetary shock persistence
eta_n = 0.8; % Nondurable consumption share
eta_d = 1-eta_n; % Durable consumption share
omega = 0.5; % Share of savers
m = 0.7; % HH LTV ratio
hab = 0.5; % Habit formation parameter
wag_rid = 0.3; % Weight on past wages in real wage friction
w_markup = 0.1; % wage markup 
epsilon_d = 11; % Differentiated good elasticity of substition for durables
epsilon_n = 11; % Differentiated good elasticity of substitution for nondurables
phi_n = 58.25; % Price adjustment cost for nondurables
phi_d = 0; % Price adjustment cost for durables
phi_pi = 1.5; % Taylor Rule coefficient on inflation
phi_y = 0; % TR coefficient on output
b_d = 1; % Production scale durables
b_n = 1; % Production scale nondurables

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model;

%% All equilibrium conditions are listed in Appendix D.1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Households

%%%%%%%%%%%%% Borrowers

%% Consumption with internal habits
eta_n*( (exp(c_n_b) - hab*exp(c_n_b(-1)))^(-sigma) - beta_b*hab*( (exp(c_n_b(+1))-hab*exp(c_n_b))^(-sigma) ) )
     = exp(lambda_b)*exp(p_n);

%% Wages (subject to ad-hoc real stickiness with markups)
exp(w) = ( nu_b*exp(h_b)^chi/exp(lambda_b) *(1+w_markup) )^(1-wag_rid) * (exp(w(-1))/exp(pi))^(wag_rid);

%% Marginal value of durable investment with second-order (investment) adjustment costs
exp(lambda_b)*exp(p_d) = exp(md_b)*( 1 -theta_c/2*( exp(c_d_b)/exp(c_d_b(-1)) -1 )^2
    -theta_c*( exp(c_d_b)/exp(c_d_b(-1)) -1 )*exp(c_d_b)/exp(c_d_b(-1)) ) 
    + beta_b*exp(md_b(+1))*theta_c*( exp(c_d_b(+1))/exp(c_d_b) -1 )*( exp(c_d_b(+1))/exp(c_d_b) )^2;

%% Value of durable holdings
exp(md_b) = eta_d*exp(d_b)^(-sigma) + beta_b*(1-delta_d)*exp(md_b(+1)) 
    + m*exp(psi)*exp(md_b);

%% Intertemporal Euler equation including borrowing constraint
exp(r) * exp(psi) = 1 - beta_b * ( ( exp(lambda_b(+1)) * exp(r) ) / ( exp(lambda_b) * exp(pi(+1)) ) );

%% Law of motion for durables 
exp(d_b) = (1-delta_d)*exp(d_b(-1)) + exp(c_d_b)*(1 - theta_c/2*( exp(c_d_b)/exp(c_d_b(-1)) -1 )^2);

%% Binding collateral constraint for borrowers
exp(r) * exp(b) = m * exp(p_d) * exp(d_b);

%% Budget for borrowers
exp(c_n_b)*exp(p_n) + exp(p_d)*exp(c_d_b) + exp(r(-1)) * exp(b(-1)) / exp(pi)
    = exp(b) + exp(w)*exp(h_b);

%%%%%%%%%%%%% Savers

%% Consumption with internal habits
eta_n*( (exp(c_n_s) - hab*exp(c_n_s(-1)))^(-sigma) - beta*hab*( exp(c_n_s(+1))-hab*exp(c_n_s) )^(-sigma) )
     = exp(lambda)*exp(p_n);

%% Wages with adjustment frictions
exp(w) = ( nu*exp(h_s)^chi/exp(lambda)*(1+w_markup) )^(1-wag_rid) * ( exp(w(-1))/exp(pi) )^(wag_rid);

%% Marginal value of durable purchases
exp(lambda)*exp(p_d) = exp(md_s)*(1-theta_c/2*(exp(c_d_s)/exp(c_d_s(-1))-1)^2
    -theta_c*(exp(c_d_s)/exp(c_d_s(-1))-1)*exp(c_d_s)/exp(c_d_s(-1))) 
    + beta*exp(md_s(+1))*theta_c*(exp(c_d_s(+1))/exp(c_d_s)-1)*(exp(c_d_s(+1))/exp(c_d_s))^2;

%% Value of durable stock
exp(md_s) = eta_d*exp(d_s)^(-sigma) + beta*(1-delta_d)*exp(md_s(+1));

%%% Law of motion for durables
exp(d_s) = (1-delta_d)*exp(d_s(-1)) + exp(c_d_s)*(1 - theta_c/2*(exp(c_d_s)/exp(c_d_s(-1)) -1 )^2);

%% Intertemporal Euler equation with no borrowing constraint
1= beta * ( ( exp(lambda(+1)) * exp(r) ) / ( exp(lambda) * exp(pi(+1)) ) );  


%%% Note: This file by default excludes saving. To add it, create another
%%%  variable s and define a new equation: omega*s + (1-omega)*exp(b) = 0. 
%%%  The s variable will not have an exp around it because unlike all other
%%%  variables in the model it will need to take negative values.

%%%%%%%%%%%%%%%%%%%%%%% Firms

%% Wage efficiency conditions
exp(w) * (1 + exp(mu)) = (1-alpha_d) * exp(mc_d) * a * b_d * exp(k_d(-1))^alpha_d * exp(h_d)^(-alpha_d);
exp(w) = (1-alpha_n) * exp(mc_n) * a * b_n * exp(k_n(-1))^alpha_n * exp(h_n)^(-alpha_n);

%% Investment efficiency conditions
exp(p_d) * (1+exp(mu)) = exp(mk_d)*(1-theta_d/2*(exp(i_d)/exp(i_d(-1))-1)^2
    -theta_d*(exp(i_d)/exp(i_d(-1))-1)*exp(i_d)/exp(i_d(-1))) 
    + beta*(exp(lambda(+1))/exp(lambda)*exp(mk_d(+1))*theta_d*(exp(i_d(+1))/exp(i_d)-1)*(exp(i_d(+1))/exp(i_d))^2);

exp(p_d) = exp(mk_n)*(1-theta_n/2*(exp(i_n)/exp(i_n(-1))-1)^2
    -theta_n*(exp(i_n)/exp(i_n(-1))-1)*exp(i_n)/exp(i_n(-1))) 
    + beta*(exp(lambda(+1))/exp(lambda)*exp(mk_n(+1))*theta_n*(exp(i_n(+1))/exp(i_n)-1)*(exp(i_n(+1))/exp(i_n))^2);

%% Capital efficiency conditions
exp(mk_d) = beta*exp(lambda(+1))/exp(lambda)*( b_d*alpha_d*exp(mc_d(+1))*a(+1)*exp(k_d)^(alpha_d-1)*exp(h_d(+1))^(1-alpha_d)
    + exp(mk_d(+1))*(1-delta) + exp(p_d(+1))*exp(pi(+1))*exp(mu(+1))*xi(+1) ) ;

exp(mk_n) = beta*exp(lambda(+1))/exp(lambda)*( b_n*alpha_n*exp(mc_n(+1))*a(+1)*exp(k_n)^(alpha_n-1)*exp(h_n(+1))^(1-alpha_n)
    + exp(mk_n(+1))*(1-delta) );

%% Binding borrowing constraint for durable producers
exp(w)*exp(h_d) + exp(p_d)*exp(i_d) = xi*exp(p_d)*exp(k_d(-1));

%% Phillips Curves
((1-epsilon_d)*exp(p_d) + epsilon_d*exp(mc_d)) - phi_d*(exp(pi_d)-1)*exp(pi_d) 
    + beta*phi_d*(exp(lambda(+1))/exp(lambda) * (exp(pi_d(+1))-1) * exp(pi_d(+1)) * exp(y_d(+1))/exp(y_d)) = 0;

((1-epsilon_n)*exp(p_n) + epsilon_n*exp(mc_n)) - phi_n*(exp(pi_n)-1)*exp(pi_n) 
    + beta*phi_n*(exp(lambda(+1))/exp(lambda) * (exp(pi_n(+1))-1) * exp(pi_n(+1)) * exp(y_n(+1))/exp(y_n)) = 0;

%% Production functions
exp(y_d) = b_d*a*exp(k_d(-1))^alpha_d*exp(h_d)^(1-alpha_d);
exp(y_n) = b_n*a*exp(k_n(-1))^alpha_n*exp(h_n)^(1-alpha_n);

%% Capital accumulation
exp(k_d) = exp(i_d)*( 1 - theta_d/2*( exp(i_d)/exp(i_d(-1)) -1 )^2 ) + (1-delta)*exp(k_d(-1));
exp(k_n) = exp(i_n)*( 1 - theta_n/2*( exp(i_n)/exp(i_n(-1)) -1 )^2 ) + (1-delta)*exp(k_n(-1));

%% Market clearing

% Goods markets
exp(c_d) + exp(i_d) + exp(i_n) + phi_d/2*(exp(pi_d)-1)^2*exp(y_d) = exp(y_d);
exp(c) + phi_n/2*(exp(pi_n)-1)^2*exp(y_n) = exp(y_n);

% Labor market
omega*exp(h_s) + (1-omega)*exp(h_b) = exp(h_d) + exp(h_n);

% Aggregates
exp(k) = exp(k_d) + exp(k_n);
exp(d) = omega*exp(d_s) + (1-omega)*exp(d_b);
exp(y) = exp(p_d)*exp(y_d) + exp(p_n)*exp(y_n);
exp(c_d) = omega*exp(c_d_s) + (1-omega)*exp(c_d_b);
exp(c) = omega*exp(c_n_s) + (1-omega)*exp(c_n_b);

%% Exogenous processes
a = a(-1)^rho*exp(e_a);
money=(money(-1))^rho_m*exp(e_m);
xi=xibar^(1-rho)* xi(-1)^rho * exp(e_xi);
nu = nu_bar*exp(e_nu);
delta = delta_bar*exp(e_delta);
sigma = sigma_bar*exp(e_sigma);

%% Taylor rule
beta * exp(r) = (beta*exp(r(-1)))^(rho) * (exp(pi)^phi_pi)^(1-rho) * money;

%% Price definitions
exp(pi_d) = exp(p_d)/exp(p_d(-1)) * exp(pi);
exp(pi_n) = exp(p_n)/exp(p_n(-1)) * exp(pi);
exp(p_d)^eta_d*exp(p_n)^eta_n=1;

%% MPK definitions (will be equal to user cost in equilibrium)
exp(mpk_d) = b_d*alpha_d*a(+1)*exp(k_d)^(alpha_d-1)*exp(h_d(+1))^(1-alpha_d);
exp(mpk_n) = b_n*alpha_n*a(+1)*exp(k_n)^(alpha_n-1)*exp(h_n(+1))^(1-alpha_n);

end;


%% Set steady state values
initval;

c_n_b    		 =-0.427943;
c_d_b    		 =-2.32963;
d_b      		 =1.5824;
h_b      		 =-0.750065;
psi      		 =-4.60517;
lambda_b 		 =0.226003;
b        		 =1.22128;
md_b     		 =0.231608;
c_n_s    		 =-0.142803;
c_d_s    		 =-1.94483;
d_s      		 =1.96719;
h_s      		 =-1.04506;
lambda   		 =-0.0689892;
md_s     		 =-0.0633844;
k        		 =2.18299;
d        		 =1.79319;
c        		 =-0.275244;
c_d      		 =-2.11883;
y        		 =0.13702;
mpk_d    		 =-3.20474;
mpk_n    		 =-3.11404;
mc_d     		 =-0.0897054;
mc_n     		 =-0.0967114;
mk_d     		 =0.0579619;
mk_n     		 =0.00560478;
y_d      		 =-0.950995;
y_n      		 =-0.275244;
h_d      		 =-2.00406;
h_n      		 =-1.28296;
k_d      		 =1.15513;
k_n      		 =1.74018;
i_d      		 =-2.35143;
i_n      		 =-1.76638;
mu       		 =-2.92337;
r        		 =0.0100503;
w        		 =0.505536;
pi       		 =0;
p_d      		 =0.00560478;
p_n      		 =-0.00140119;
pi_d     		 =0;
pi_n     		 =0;
a        		 =1;
money    		 =1;
xi       		 =0.1;
nu       		 =4;
delta    		 =0.03;
sigma    		 =1;


end;


% This is the default steady state solver if you don't use homotopy
steady(solve_algo=1,maxit=100000);

%%%% If you want to use homotopy, use this block =======
%
%homotopy_setup;
%
% Change parameters here using the following formatting:
%phi_pi, 1.5; 
%
%end;
%steady(homotopy_mode=1,homotopy_steps=100,solve_algo=1);
%=======================================================


check;

% Specify which shocks to include
shocks;
var e_m = 0.01;
var e_a = 0.05;
%var e_nu = 1;
%%var e_g = 0;
%var e_sigma = 1;
var e_xi = 0.1;
%var e_delta = 0.1;


end;

%% Set Dynare seed. This won't matter for linearized IRFs, but 
%%  it can potentially matter for simulations
set_dynare_seed(123);

% Optional: Run model diagnostics if something isn't working
%model_diagnostics(M_,options_,oo_)

%% Stochastic simulation
stoch_simul(periods=50000, pruning, irf=17, order=1, nograph) c k_d k_n k y pi p_d i_d i_n mpk_d mpk_n mu;

%%%%% Save output for plotting

% Extract variable names (format is variable_shock, e.g. k_d_e_m)
var_names = fieldnames(oo_.irfs);
shock_irfs = zeros(length(var_names),17);

% Loop through the structure to get all IRFs in one matrix
for i=1:length(var_names)
shock_irfs(i,:) = oo_.irfs.(var_names{i});
end;

%% Combine IRFs with names of variable and shock
shock_irfs = [var_names string(shock_irfs)];

%% Save output as .csv file to be plotted 
writematrix(shock_irfs,"irfs\shock_irfs_baseline.csv")



